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ABSTRACT 

A graphical test bed in which the results of a simulation 
experiment can be reported and analyzed is described. The 
test bed is based on the regression adjusted graphics and 
estimation (RAGE) methodology developed by Heidelberger and 
Lewis [Ref. 1] for regenerative simulations. From the graphics 
and the associated numerics the experimenter can summarize 
and see simultaneously relative properties, such as bias, 
normality and standard deviation, of several estimators of a 
characteristic of a population for up to eight sample sizes. 

The graphics is supported on a line printer to make it and 
the program portable. 
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I. INTRODUCTION 



During the course of a simulation experiment the experi- 
menter is frequently confronted with the problem of estimating 
some function or parameter from the generated data. He will 
usually have a choice of one or more estimators, parametric 
or nonparametric , that can be used to estimate the function 
or parameter. These estimators will vary in their statistical 
properties. Some will estimate with low variance but contain 
bias, others will have just the opposite characteristics and 
still others will contain varying degrees of both. How these 
different estimators compare in relation to each other as 
well as over varying sample sizes are of concern to the experi- 
menter. Comparison of estimators are important in that they 
help the experimenter to choose the most appropriate estimator 
for his particular analysis. Comparisons may also be helpful 
in deciding whether a new estimator is of any use. The appro- 
priate estimator to use is a function of the bias, variance 
and nonnormality in the estimate for available run length 
(sample size) . 

When the simulated data from which the properties of the 
estimators is to be examined is independent and identically 
distributed (iid) it can be sectioned into M independent 
blocks of specified sample size (N) . From each of these 
blocks the estimator can be applied to all subsamples of 
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size , thereby creating samples of estimates of varying 

size for each estimator. As a first cut at the properties 

of the estimators, box plots of these samples of estimates 

are drawn for subsample sizes n (see, e.g. Figure 3a). 

k 

However, something is usually known about the properties 

of the estimators and these properties can be used to enhance 

the graph of box plots versus subsample size, n^,. Thus it 

is often known that bias in the estimator is a polynomial 

in 1/n, and the variance of the estimator is a polynomial in 

1 , 1 , 1 , ... Also the estimators are usually asymp- 

n n l * 5 n^ 

totically normally distributed. Consequently the graphs are 
enhanced by an estimated unbiased regression line for the 
mean of the estimates of the parameter as well as a horizontal 
line for the regression-estimated expected value of the 
parameter (asymptote) . A regression is also performed for 
the variance structure of the estimators and the sample 
characteristics (kurtosis and skewness) are given under each 
box plot to assist in judging the convergence to normality. 

Heidelberger and Lewis [Ref. 1] used this methodology in 
the context of regenerative simulation to develop a procedure 
and protocol to produce unbiased and normally distributed 
regression adjusted regenerative estimates (RAGE), ;vith a 
given precision or for total available computing time. How- 
ever, the graphics associated with this procedure required a 
Tektronic 4013 terminal using an APL graphics package making 



12 



its use limited to those organizations possessing such a 
terminal. The protocol also requires a relatively experi- 
enced user to properly interpret the graphical and statisti- 
cal output and to apply the sequential steps of the protocol 

The purpose of this thesis is to apply the methodology 
developed by Heidelberger and Lewis [Ref. 1] to comparing 
the bias and variance of estimators of a parameter from simu 
lated output. This objective also includes supporting the 
graphics associated with the methodology on a line printer 
to make its use more accessible. 

Details of the graphics and methodology are given in 
Chapter II, entitled USAGE. Details of the arguments of the 
RAGE program are given in Chapter III and a sample program 
using the RAGE graphics is given in Chapter IV. All the 
user has to supply is the input data and FORTRAN functions 
to compute the estimators for any given subset of data. 

A comparative study of estimators of the gamma shape 
parameter using the graphical test bed RAGE is given in 
Chapter V. This comparative study illustrates the different 
properties of estimators that can be seen with regression 
adjusted graphics as sample size increases. 
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II. USAGE 



The regression adjusted graphics and estimating (RAGE) 
program is designed to evaluate distributional properties of 
one or more estimators as a function of sample size from a 
stationary series of raw data. These estimators are program- 
med as functions by the user and declared external in the 
main program. The graphical and numerical results that 
accompany each estimator are varied and will depend to some 
extent on the options desired by the user. The purpose of 
the following discussion is to explain the graphical and 
numerical output as well as how they were obtained and some 
of the options available to the user. 

The general procedure begins with a given iid data base 
or one that has been simulated and an unknown parameter or 
population characteristic, 0, that is to be estimated by an 

A 

estimator, 0, from the data. The experimenter must first 

decide on the number of independent blocks or sections (M) 

of size N the data is to be sectioned into. These M sections 

or blocks will henceforth be referred to as replications; 

the terminology will seem more appropriate as the discussion 

* 

continues. The number of replications (M) can alternatively 
be determined indirectly by the sample size (N) desired in 
each section. In either case MxN must not exceed the total 
allowable sample size. This sectioning of the data into 
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M replications of sample size N is represented in Figure 1 
where the entire block represents the MxN samples from the 
total sample population. 



N 
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m. 



Tu- 
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i i=l, 



M 



N 



m 



M 



Figure 1: M replications, each of sample size N 

Within each replication the experimenter must decide on 

/\ 

the subsample sizes that the estimator, 0, is to use in 
estimating the parameter. The number of different subsample 
sizes, K, that can be used to resection the same data can vary 
from one to eight and the size of each subsample can vary 
from one to the size of the sample, N. This repeated sub- 
sectioning of the same data is represented pictorially by 
Figure 2 where [x] refers to the highest integer less than 
or equal to x and r^, k=l , . . . , K refers to the number of 
estimates possible from a sample of size N using subsample 
sizes nj., k=l, ..., K to estimate the parameter. Each box 
or rectangle in Figure 2 represents the same MxN data set 
with the same M replications of size N subsectioned into 
different sections of size nj,. Note that the n^'s divide 
N evenly in Figure 2. This is not necessary but certainly 
desirable since all available information will be utilized. 
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Figure 2: Subsectioning MxN samples into subsamples 

of size n k , k=l , . . . , K 

If each subsample size (n k , k=l, . K) produces a total 
of r k , k=l , • • • , K estimates per sample of size N then there 
will be a population of estimates from the M replications 
that total Mxr k , k=l, K for each subsample. For each 
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population of estimates from each subsample (n^., k=l, K) 

box plots are plotted so that the distribution of the esti- 
mates for the varying subsample sizes can be seen graphically. 

In addition to these box plots of the estimates, a 
regression adjusted estimate (RARE) of the parameter is 
computed and plotted over the boxplots. In addition to the 
RARE estimate, unbiased RARE estimates out to a specified 
number of terms determined by the user are computed and plotted 
on the graph (The mathematical derivation of the RARE and 
unbiased RARE estimates are given later in this section) . 

The unbiased RARE estimates are computed as a function of 
subsample size and plotted in conjunction with the RARE 
estimate so that a comparison of the two values can be made 
graphically and the evolution to an unbiased state can be 
seen as the subsample size increases. 

This idea can be made clearer through the use of an example 
and its associated graphics and numerics, Figure 4a. Here an 
independent and identically distributed sample of 50,000 
gamma distributed random numbers (shape parameter, K=5.0) 
has been generated and sectioned into M=100 sections (replica- 
tions) with each replication containing N=500 randon numbers 
(sample size). This information is printed at the top of 
the graph for future identification. The experimenter chose 
eight different subsample sizes at which to examine the 
parameter (33, 50, 71, 125, 166, 250 and 500). These 
subsample sizes (n^'s) are printed on the horizontal scale 
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of the graph. The population of estimates of the parameter 
at each subsample size are then plotted as vertical box plots 
above these numbers. When the population of estimates deter- 
mined from a subsample size is less than nine, box plots 
are omitted and the values of these estimates are simply 
plotted with the sumbol "0" on the graph. There are no 
estimates with populations totaling less than nine in Figure 
4a but the graphical output for such cases can be seen in 
Figure 3a where subsample sizes, 125, 166, 250 and 500 total 
less than nine estimates. 

The horizontal scale of the graph is scaled roughly to 
the minimum and maximum subsample sizes with the box plot of 
each subsample's estimates spaced proportionally. However, 
when subsample size differences are relatively small so that 
overlap of the box plots would occur, a minimum of six spaces 
between centers of box plots is left to prevent the graphics 
from becoming confused. 

The vertical scale of the graph is determined by the 
minimum and maximum values from the population of estimates 
determined from the smallest subsample size (n^) . It is 
assumed here that subsequent subsample sizes (n ? , n-, ...) 
will be large in relation to the smallest subsample size 
(n^) so that estimates produced from these will be within 
the range of that determined from the smallest subsample 
size. This assumption appears tenable if the range of esti- 
mates is compared in Figure 4a for increasing subsample sizes. 
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The vertical scale of Figure 4a is not scaled to the 
minimum and maximum value of the estimates determined from 
the smallest subsample size (n^=33) . This is because the 
RAGE program can evaluate up to three different estimators 
at one time and txvo alternate estimators (Figures 4b 5 4c) 
were also used to estimate the parameter. An option was 
entered to scale all three estimators to the same scale and 
it chose the minimum and maximum values from all the esti- 
mates produced by all three estimators using the smallest 
subsample size (n-^33). By comparing Figures 4a, 4b and 4c 
it can be seen that the maximum value was taken from the 
first estimator. Figure 4a, and the minimum value was taken 
from the third estimator, Figure 4c. Other scaling options 
will be discussed in detail in the following chapter. 

Below the graph statistical results as they pertain to 
the population of estimates for each subsample size are 
printed. These statistics include the sample mean, standard 
deviation, standard deviation of the mean, skewness and kur- 
tosis. Thus, for example, in Figure 4a at subsample size 
n-^ = 33 there is a population of [N/n-jJxM = r^xM = [500/33]xl00 
150 estimates for which a sample mean, standard deviation, 
standard deviation of the mean, skewness and kurtosis is 
computed. The subsample sizes these statistics refer to are 
labeled across the top. 

As mentioned earlier, in addition to box plots and associ- 
ated statistics of the estimates, RARE and unbiased RARE 
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estimates of the parameter are computed and plotted on the 
graph. The degree of regression, D, for the computations of 
the unbiased RARE estimates is determined by the user and 
printed at the top of the graph for future reference. The 
RARE estimate is the regression coefficient 3 q and is the 
estimated asymptotic expected value of the parameter being 
estimated. The unbiased RARE estimates are determined by 
the regression coefficients, 3-^, $ 2 > •••> 3^, the subsample 
size, n, and the RARE estimate, 3q, using the following 
formula : 



A set of regression coefficients (3q, 3^, ...» 3^) is 
determined from the mean of the population of estimates 
(dependent variable) at each subsample size for each replica 
tion. The M regressions that determine the M sets of coeffi 
cients (Bq, 3-j, > •••> 3^) are represented by the following 

formula : 





m" ( n d = 6 0 




1 




n 



1 



1 



m 3(n k ) - 3 0 + S]_ _ 





m = 1 



• • • 9 



M 
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where § (n. ) 
m K 




k = 1, . . . , K 
m = 1 , . . . , M 



The regression coefficients, to determine the RARE and unbiased 
RARE estimates, are the average values determined from these 
M regressions applied to the M replications. The mean, 
variance and standard deviation of these coefficients are 
printed below the graph below the subsample site statistics. 

As mentioned earlier, the first coefficient, $q, is the 
estimated asymptotic expected value of the parameter being 
estimated (RARE). Its mean from the M regressions is printed 
as a horizontal dashed line across the graph. As the sub- 
sample size increases the bias in the unbiased RARE estimates 
becomes smaller and the unbiased RARE estimates (dotted line) 
will asymptotically approach the RARE estimate (dashed line) . 

In Figure 4a the estimated expected value or asymptote 

8 q = 4.98 and the bias coefficients for the three degree (D) 

regression are: 3]_ = 16.50 , 32 = "182.21, 83 = 3659.04. 

The variance regression coefficients are determined from 
the variances of the population of estimates for each sub- 
sample size (dependent variable) and evaluating the single 
regression 
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VAR(0 (n, ) ) = a x _L_ + * J_ + 

1 n-. L 1.5 

1 n. 



a 



d+1 T d + 2)/2 
n. 
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m 



M r k 

S A 

= 1 i = l m'’i^ n k^ 



k=l , . . . , K 



r k M 



k=l , .... K 



The four (D+1) variance regression coefficients for Figure 4a 
are; a-^ = 58 . 37 , a 7 = 174.72 , a-j = -1500.13 , and = 5446.25. 
The variance regression coefficients derived here can be used 
to obtain an approximate variance for a specified subsample 
size by substituting the subsample size for n^ in any one 
of the above regression equations along with the variance 
coefficients (a^ , a ? , ...) and completing the computation. 

Below the variance regression coefficients the minimum 
and maximum value of the vertical scale is printed. The 
purpose of this information is to allow the user to scale 
other estimators of the same parameter to the same scale. 

It can be used when the number of estimators exceeds three 
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(the maximum number of estimators RAGE can use in a single 
run) or when subsequent runs for different conditions are 
desired to be of the same scale for graphical comparisons. 

A final identifying feature of the graphics is a section 
at the very bottom of the graphical output which allows the 
user to identify the estimator that was used in estimating 
the parameter. In Figure 4a this estimator was the moment 
estimator, i.e. the reciprocal of the estimated coefficient 
of variation squared. 
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III. DESCRIPTION OF RAGE ARGUMENTS 



The RAGE program has the calling sequence: call 

RAGE (X ,N ,M ,NE , L , D , RG , SEI , SVS , YMIN , YMAX , NEST , EST1 , TTL1 ,EST2 , 
TTL2 ,EST3 ,TTL3) . The arguments are relatively easy to 
understand but in several cases are related and must meet 
certain restrictions. For this reason a detailed descrip- 
tion of each argument is given, as well as any restrictions 
or limitations that must be met within or between arguments. 
Following this detailed discussion an example is given in 
Chapter IV on how to set up the main program and the estima- 
tor functions. 

A. DATA ARRAY ... X 

The X argument is a single precision Real*4 array con- 
taining the data from which the parameter is to be estimated. 
The only restriction is that the size of the X array must 
not exceed 50,000 which is the maximum sample size that RAGE 
can store. 

B. SAMPLE SIZE. . . N 

The number of data elements per section or sample size 
is designated by the integer N. Depending on the minimum 
sample size required by the estimator, N can vary from 1 to 
50,000. The sample size, N, times the number of sections, 

M, cannot exceed 50,000 (the maximum number of data elements 
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stored in the array X). If M times N exceeds 50,000 an 
error message will be printed and execution will terminate. 

If M times N exceeds the number of data elements stored by 
the user in the array X then no error message will be printed 
and the program will execute giving unpredictable results or 
an abnormal ending will occur. 

C. NUMBER OF REPLI CATIONS ... M 

The number of sections or replications the array X is 
sectioned into is determined by the integer value assigned to 
M. The number assigned to M is also the number of regressions 
that will be run to determine the regression coefficients 
Bq> B , ..8^. The number of regressions is limited to one 
hundred and therefore M is limited to a maximum value of one 
hundred. An error message will be printed and execution will 
terminate if M is greater than one hundred or less than one. 

If M is equal to one the variance and standard deviation of 
the regression coefficients for the unbiased estimate will 
not be printed since they will not exist. 

D. SUBSAMPLE SIZE ARRAY... NE 

The argument NE is an integer array of size eight con- 
taining the subsample sizes (n^, k=l, ..., K) the estimator 
is to use in estimating the parameter. The values in the 
array NE should be ordered from smallest to largest or a 
warning message will be issued to indicate the vertical 
scale of the graph may not be large enough to include all 
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estimates. RAGE will use the population of estimates deter- 
mined from the subsample size stored in NE(1) to determine 
the vertical scale of the graph. If NE(1) is not the smallest 
value the range of the graphics may not be large enough to 
include outliers determined from subsequent smaller subsample 
sizes . 

Due to storage limitations a maximum of 12,500 estimates 
is permissable. Since NE(1) will produce the most estimates, 
the user can determine the total number of estimates using the 
formula Mx[N/NE(l)]. If this value exceeds 12,500 an error 
message will result stating this limitation has been exceeded 
and execution will terminate. Corrective action can be taken 
by increasing NE(1) or decreasing M or N. 

E. NUMBER OF SUBSAMPLE SIZES (BOX PLOTS)... L 

The user can select the number of subsample sizes to 
be used from the array NE(8) by assigning an integer value 
to the argument L. The argument, L, can take on values 
between one and eight. RAGE will then select the first L 
subsample sizes from NE(8), compute the population of esti- 
mates of the parameter with each, graph the box plots and 
print the associated statistics. If L is not an integer 
value between one and eight an error message will result and 
execution will terminate. 

If L = 1 no regressions are possible and any regression 
output is omitted. The maximum degree of regression allowed 
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by RAGE, regardless of the value assigned to the argument D 
(degree of regression), is L-l. If D exceeds this value RAGE 
will set the degree of regression to L-l. 

F. DEGREE OF REGRESSION. . .D 

The degree of regression, D, is an integer value less 
than or equal to three. The degree of regression refers to 
the number of bias coefficients, 8 ]_, 67 , ..$ 3 , that will be 
used in computing the unbiased estimate out to D terms. 

In addition to the restrictions on D mentioned earlier, 

D will also be reduced for the regression that determines 
the variance coefficients, a^, a 7 , •• a c [+i> when the number of 
subsamples sites with a variance is less than D+l. 

G. REDUCED GRAPHICS. . . RG 

The reduced graphics argument, RG , allows the user to 
reduce the range of the vertical scale of the graph by elimi- 
nating outliers and only printing the number of outliers 
that exist. The vertical scale is reduced to the upper quar- 
tile plus 1.5 times the interquartile distance value and the 
lower quartile minus 1.5 times the interquartile distance 
value of the estimates from the smallest subsample site, 
NE(1). The vertical scale then becomes these values and the 
outliers are merely counted and printed on the graph at the 
ends of each box plot's upper (lower) quartile + (-) 1.5 
times interquartile distance position. If there are no out- 
liers for the NE(1) estimates, either or both extremes, the 
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vertical scale is further reduced to the first estimate 
value within these limits. RG = 1 will cause reduced graphics 
to be produced and RG = 0 will cause full scale graphics. 

H. SCALING ESTIMATORS INDIVIDUALLY ... SEI 

The argument SEI allows the user to scale each estimator's 
vertical scale for the graphics separately or to the same 
scale. SEI = 0 causes all estimators' vertical scale to have 
the same range of values. The range of the vertical scale 
and whether outliers are printed or counted will depend on the 
graphics desired (RG = 0 or 1). Figures 4a, 4b and 4c are the 
results of three estimators scaled to the same scale (SEI = 0) 
with full scale graphics (RG = 0). Figures 5a, 5b and 5c are 
the results of three estimators scaled to the same scale 
(SEI = 0) but with reduced graphics (RG = 1) . 

If SEI = 1 each estimator's vertical scale will be scaled 
individually. The range of the vertical scale and whether 
outliers are printed or counted will depend on the graphics 
desired (RG = 0 or 1) . 

I. SETTING THE VERTICAL SCALE... SVS 

To set the vertical scale an integer value of 1 is assigned 
to the argument SVS. SVS = 1 will cause the SEI argument to 
be ignored and the vertical scale of all the estimators' 
graphics will be set to those values assigned to YMIN and 
YMAX. When the argument SVS is in effect (SVS = 1) the 
argument RG is also in effect as far as whether outliers 
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are printed or counted (RG = 0 or 1) but the determination 
of the vertical scale by RG is overridden by SVS and is set 
to YMIN and YMAX. SVS = 0 causes RAGE to ignore values 
assigned to YMIN and YMAX and scaling is determined by the 
values assigned to the arguments RG and SEI. Figure 4d is 
an example of SVS = 1 and RG = 0 with the vertical scale 
set to that of Figures 4a, 4b, and 4c (YMIN=-1 . 2812 , 

YMAX=13 . 9984) . Figure 5d is an example of SVS = 1 and 
RG = 1 with the vertical scale set to that of Figures 5a, 

5b and 5c (YMIN=1 . 0394 , YMAX=9 . 1423) . 

J. MINIMUM VALUE OF THE VERTICAL SCALE... YMIN 

The value YMIN is a real number assigned by the user as 
the lower limit of the vertical scale. If YMIN is assigned 
a value that is so large that it falls within the range of 
values that determine the body of any box plot then an error 
will occur and an abnormal ending will result. 

K. MAXIMUM VALUE OF THE VERTICAL SCALE... YMAX 

The value YMAX is a real number assigned by the user as 
the upper limit of the vertical scale. If YMAX is assigned 
a value that is so small that it falls within the range of 
values that determine the body of any box plot then an error 
will occur and an abnormal ending will result. 

L. NUMBER OF ESTIMATORS ... NEST 

The argument NEST is an integer, 1, 2 or 3, which indicates 
the number of estimators that will be used to estimate the 
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parameter. An error message will be issued if NEST is not 
one of these values and execution will terminate. 

M. ESTIMATOR FUNCTIONS ... EST1 , EST2 , EST3 

These arguments are the function names assigned to each 
of the estimators that are to be used in estimating the para 
meter. Each function must be written with a two argument 
calling sequence (i.e. call EST1(X,N)) where the first 
argument, X, must be the data array and the second argument, 

N, must be the number of data elements. 

If there are less than three estimators then dummy argu- 
ments must be inserted in the calling sequence where subse- 
quent estimators would normally appear (Suggestion .. .use a 
function of a previous estimator, i.e. EST1, as the dummy. 
RAGE will not provide duplicate output for these repeated 
functions unless NEST is increased to include them.). Each 
function must then be declared external in the main program. 
For an illustration of how the functions are written and 
declared see the example in Chapter IV. 

N. DESCRIPTION OF ESTIMATORS .. .TTL1 ,TTL2 ,TTL3 

The titles associated with each estimator can be passed 
as an argument through the call program (RAGE) provided 
there is a minimum of 120 characters between apostrophes 

(i.e. call RAGE(..., EST1,' 120 characters ' ,EST2,...)) 

An alternate and possibly easier method is to declare each 
title (TTL1, TTL2, TTL3) as a Real*8(15) array and assign 
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the characters that describe each estimator as data to each 
title. For an example of this method see the example in 
Chapter IV. 

The user should note that the title associated with each 
estimator must immediately follow the estimator argument. 

As in the case of estimators, dummy titles must be inserted 
when the number of estimators requiring titles is less than 
three. Finally, it should be noted that the program will 
not execute properly unless titles of the correct format are 
included. 
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IV. EXAMPLE OF A PROGRAM USING RAGE 



In the following example suppose a user has a set of iid 
data from which he wishes to estimate the population mean by 
using two estimators, the sample mean and the sample median. 
Assume that these estimators are given the abbreviated func- 
tion names MN and MED and that the associated titles for 
each estimator are "The sample mean as an estimator of the 
population mean" and "The sample median as an estimator of 
the population mean - ." If the subsample sizes at which it is 
desired to examine the properties of these estimators are 
33, 50, 71, 125, 166, 200, 250 and 500 the FORTRAN program 
may look as formatted on the following page. 
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REAL *4 X (50000) ,YMIN,YMAX 
REAL*8 T1 (15) ,T2 (15) 

INTEGER NEST ,D ,RG ,SEI ,SVS ,M,N ,L ,NE (8) 

DATA NE/ 33,50,71,125,166,200,250 500/ 

INSERT 8 CHARACTERS BETWEEN EACH SET OF APOSTROPHES UNTIL 
TITLE IS COMPLETED AND FILL REMAINDER WITH BLANKS. 

DATA Tl/'THE SAMP ' , ' LE MEAN ' , ’AS AN ES’,'TIMATOR ', 

+ ’ OF THE P ' , ' OPULATIO ’ , ' N MEAN', 8*' '/ 

DATA T2/ ' THE SAMP ' , ' LE MEDIA' , ' N AS AN ’ , ' ESTIMATO ' , 

+ 'R OF THE' , 'POPULATI ’ , 'ON MEAN ' , 8 * ' '/ 

EXTERNAL MN ,MED 



CALL RAGE (X ,N ,M ,NE , L ,D , RG , SEI ,SVS , YMIN , YMAX , 
+NEST ,MN ,T1 ,MED,T2 ,MN,T1) 



STOP 



REAL FUNCTION MN(X,N) 

REAL X (N) 

C AVERAGE OF N VALUES IN X 

MN= 0 . 0 
DO 10 1 = 1, N 

MN = MN +X(I) 

10 CONTINUE 

MN = MN/ FLOAT (N) 

RETURN 

END 

REAL FUNCTION MED(X ,N) 

REAL X (N) 

C MEDIAN OF N VALUES IN X 

C "PCTILE" IS A FUNCTION THAT RETURNS THE 50 PERCENTILE 

C NUMBER FROM THE DATA ARRAY X OF SIZE N. 

MED = PCTILE (X,N, . 5) 

RETURN 

END 

In this example since only two estimators were used a dummy 
estimator and title had to be inserted as the arguments for 
the third case. Here the first estimator and its title were 
simply repeated. 
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V. A COMPARATIVE STUDY OF ESTIMATORS FOR THE GAMMA SHAPE 



PARAMETER 



To demonstrate the use of the Regression Adjusted Graphics 
and Estimating (RAGE) package, two arrays of gamma distri- 
buted random variables with shape K=5.0 and K=0.2 were 
generated using the Naval Postgraduate School's random number 
generator- -LLRANDOMI I [Ref. 3]. Because these random numbers 
are independent and were drawn from the same distribution, 
they can be sectioned into independent and identically 
distributed blocks of size N. RAGE was specifically designed 
to evaluate estimators of parameters from just such iid 
populations. These two samples then should demonstrate the 
capabilities of the RAGE package to evaluate the bias, vari- 
ance and distribution of estimators of a specified parameter 
from the populations. 

In the following examples it is the shape parameter, K, 
of the gamma distribution that is being estimated. There are 
two primary estimators used in these examples. The first is 
the moment estimator (reciprocal of the estimated coefficient 

1 y 2 

of variation squared), K = = £ . The second 

(C(X)) 2 VAR (X) 

estimator is the maximum likelihood estimate (MLE) [Refs. 2,4] 
which is determined by a bisection search for the shape para- 
meter, K, until the following equation is satisfied: 
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n 

5" Ln X. 

Ln(k) - Digamma(K) = Ln X - iz]: i = l, . .., n. 

n 

The X^, i=l, 2, n are the random numbers from a specified 

block of size n. The other two estimators investigated are 
jacknifed versions of the initial estimators. The jacknife 
routine [Ref. 5] with each of these estimators is a four fold 
jacknife which has the form 

a- % a. ^ 

where 0*^ = 4(0 a -j^) - 3(0j), 0 a n is the estimate from the 

a, 

entire sample (n) and 0j is the estimate with the j t ^ group 
missing. Each group is a set of one fourth of the data 
points . 

Each of the four estimators is examined for varying sub- 
sample sizes and replications to determine how accurately 
each estimator estimates the shape parameter of the gamma 
distribution as the subsample size increases. The shape 
parameter of the gamma distribution was chosen as an example 
because the behavior of the many estimators which have been 
proposed for K is fairly well known. Roughly speaking, it is 
known that the moment estimator is quite efficient for K>2 
but quite inefficient for K<1 even for large sample sizes. 

Bias and normality of these estimators is not documented. 

Again the jacknife is known to work well with moment estimators 
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even in small samples but no such results are known for the 
jacknifed MLE estimator except that it is asymptotically 
efficient . 



For this investigation there are basically three examples 
representing different simulation scenarios for evaluating 
each of two cases, K=5.0 and K=0.2. They are: 

(a) Small subsample sites (33 to 500) , small 
replications (M=2) ; 

(b) Small subsample sizes (33 to 500) , large 
replications (M=100) ; 

(c) Large subsample sizes (100 to 2000) , large 
replications (M=25). 

The first can be looked on as a pilot study to determine 
regions of interest for the two longer scale studies. How- 
ever, the procedures are not fully sequential, unlike the 
usage of the RAGE methodology proposed by Heidelberger and 
Lewis [Ref. 1] for regenerative simulations. 

With many of these examples another set of reduced graphics 
is included which reduces the vertical scale and provides 
graphics which are more suitable for examining the bias 
structure of the estimators. The reduced graphics truncates 
the outliers in the box plots giving only their number and 
not their position. 
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A. EXAMINING THE GAMMA DISTRIBUTION FOR THE SHAPE PARAMETER 
(K= 5 . 0 ) 

1 . Case 1, Example 1: Small Subsample Sizes and 

Small Replications (K=5.Q)‘ 

Figures 3a, 3b, 3c and 3d show the graphical and 
statistical results from estimating the shape parameter with 
the four different estimators. These graphs all have the 
same vertical scale so that meaningful comparisons can be 
made of the box plots and of the bias of each estimator. 

One of the characteristics of the MLE estimator is 
that its estimates will asymptotically have the smallest vari- 
ance of all estimators. For small subsample sizes (such as 
the case of Figures 3a § 3b) the asymptotically smaller vari- 
ance of the MLE estimator is not apparent in the box plots. 

The body of the box plot for the MLE estimator for n^ = 33 
is smaller than the same box plot for the moment estimator. 
However, the MLE has more and larger outliers, indicating 
that the MLE is more nonnormal at this sample size than the 
moment estimator. If, however, the estimated standard devia- 
tion of the various subsample sizes are compared for each of 
the two primary estimators (Table 1) it can be seen that the 
MLE estimator consistently produces a smaller variance than 
does the moment estimator. From Figures 3a § 3b the bias in 
the two estimators appears comparable, starting high but 
becoming negligible by subsample size 250. 
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Table 1: Standard deviations of estimates of the shape 

parameter using the moment and MLE estimators 
for small sample sizes (replications = 2) 

Standard Deviation 



SAMPLE SIZE 


MOMENT ESTIMATOR 


MLE ESTIMATOR 


33 


1.517 


1.306 


50 


.9942 


.7151 


71 


. 8829 


. 6051 


100 


.6312 


.4612 


125 


. 7725 


.4823 


166 


.8310 


.5114 


250 


.4660 


.2772 


500 


. 5589 


. 3087 



Jacknifing is a technique that can be used on iid sam- 
ples to reduce the bias in the estimates caused by the estimator 
[Ref. 5]. An estimator with a complete jacknife is asymptoti- 
cally normal with the same variance as the unjacknifed estimator, 
but small sample properties are difficult to obtain analyti- 
cally. Figures 3c 6 3d are the jacknifed estimators and 
again because of the small number of estimates it is not 
apparent that the jacknife is contributing to nonnormality 
or increasing variance when the box plots are compared with 
those of Figures 3a 5 3b. The increased normality, however, 
can be seen in the jacknifed moment estimator by comparing 
box plots (Figures 3a 6 3c) and the coefficients of skewness 
(Table 2) with the nonjacknifed form. The change in normality 
or symmetry of the estimates from the jacknifed MLE estimator 
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is not discernable in a comparison of the box plots (Figures 
3b § 3d) or the coefficients of skewness (Table 3). 



Table 2: Skewness coefficients of the estimates of 

the shape parameter using the moment and 
jacknifed moment estimators (replications = 2) 

Skewness Coefficients 

SAMPLE SIZE MOMENT ESTIMATOR JACKNIFED FORM 



33 


.553 


. 327 


50 


. 076 


- .348 


71 


.339 


.253 


100 


.499 


.157 


125 


1.573 


.724 


166 


1.018 


.501 


250 


-.150 


- .235 


500 


N/A 


N/A 



Table 3: Skewness coefficients of the estimates of the 

shape parameter using the MLE and jacknifed MLE 
estimators (replications = 2) 

Skewness Coefficients 

SAMPLE SIZE MLE ESTIMATOR JACKNIFED FORM 



33 


1 .198 


.050 


50 


.510 


. 616 


71 


.193 


- .163 


100 


.773 


. 643 


125 


1.451 


1.216 


166 


1 .101 


1.229 


250 


- . 030 


. 081 


500 


N/A 


N/A 
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The bias in each estimator is depicted by the space 
between the dotted line and dashed line of each of the graphs. 
When the unbiased estimate (dotted line) approaches the 
estimated expected value (dashed line) then the bias becomes 
insignificant. Each estimator can be expected to vary in 
the amount of bias it has in its estimates and this bias can 
be expected to diminish as sample size increases. Comparing 
Figures 3a 5 3b it appears that the MLE estimator has about 
the same bias as the moment estimator but this is hard to 
determine at this point since the bias coefficients were 
determined from a small number of estimates. In the esti- 
mates from the jacknifed estimators one would expect the bias 
to be reduced because jacknifing is specifically designed to 
remove the 8 ^ ( 1 /n) term in the bias of the estimates but 
this is not obvious from the graphical comparisons (Figures 
3a § 3c and Figures 3b 5 3d). For the case of the MLE esti- 
mator the coefficient of the 1 /n term, 3 -^, has actually 
increased from 32.65 to 75.33 after jacknifing. The standard 
deviation of these bias terms are so large (19.14 and 16.23 
respectively), though, that they suggest more estimates 
(increased replications) are necessary to better determine 
the bias. 

Determining how much more bias one estimator has than 
another can also be accomplished by comparing the bias 
regression coefficients 3j_ , 82 ? etc. The larger these coef- 
ficients are the greater is the bias of that particular 
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estimator. Comparing the coefficients for the moment and 

MLE estimators of Figures 3a § 3b show that the MLE estimator 

has a larger bias. For a subsample size of 500 and a three 

degree regression the bias for the moment estimator is 

7.40128 183.099 1947.86 

+ + = .0155 

500 500 2 5 0 0 3 

and the bias for the MLE estimator is 

32.6486 1998.68 47308.0 

- + = .0577 

500 5 00 2 500 3 

These comparisons seem to indicate that the MLE estimator 
has more bias. This is not necessarily true, though, for 
reasons mentioned earlier. Increasing the number of replica- 
tions will alleviate much of this uncertainty since more 
precision in the estimates of the bias coefficients will 
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ESTIMATOR: FOUR FOLD J ACKNl FED MOMENT ESTIMAIUR OF THE SHAPE PARAMETER OF THE GAMMA DISTRIBUTION 

Figure 3c 
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2 . Case 1, Example 2: Small Subsample Sizes and Large 

Replications (K=5.0) 

It was pointed out in Example 1 of this case that 
the small number of replications left some doubt as to the 
accuracy of the estimates of the bias coefficients as well 
as the normality of the estimates. Figures 4a to 4d (cor- 
responding to Figures 3a through 3d) are for increased 
replications (M=100) . This has the effect of increasing the 
population of estimates for each subsample size (n^., k=l, 

K) and the number of estimates of the bias coefficients 
(3 1 , $ 2 , •••) but holds the subsample size at which the esti- 
mators are examined small as before. On first observing the 
graphs, it becomes apparent from the box plots of the two 
smallest subsample 'sizes (n-^ = 33 and n 2 = 50) that the estimates 
are nonnormally distributed (Figures 4a § 4b) . The jacknifed 
estimators remove some of this nonnormality (Figures 4c § 4c) , 
especially for the jacknifed moment estimator. In effect, 
outliers appear now at both high and low values but consider- 
able skewness still remains. 

Here it is apparent by the statistical output that the 
MLE estimator has significantly less variance in its popula- 
tion of estimates than does the moment estimator (Table 4) . 
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Table 4: Standard deviations of the estimates of the 

shape parameter using the moment and MLE 
estimators for varying sample sizes 



(replications = 


100) 






Standard 


Deviation 


SAMPLE SIZE MOMENT 


ESTIMATOR 


MLE ESTIMATOR 


33 1 


.486 


1 . 396 


50 1 


.175 


1.075 


71 


.9818 


. 8549 


100 


.8008 


.6636 


125 


.7242 


. 5918 


166 


.6499 


. 5145 


250 


. 5082 


.4036 


500 


.3417 


. 2730 



The increased symmetry and normality by the jacknifed moment 
estimator is graphically visible if Figures 4a 5 4c are com- 
pared. The change in symmetry or normality of the estimates 
of the jacknifed MLE estimator was not discernable in the 
graphics (Figures 3b § 3d) or skewness coefficients (Table 3) 
for Example 1. However, for increased replications, (M=100) 
if the box plots of Figures 4b 5 4d are compared it appears 
that the population of MLE estimates are more normally 
distributed (outliers at both ends) after jacknifing. The 
movement toward symmetry is more obvious by comparing the 
skewness coefficients of the MLE estimates and the jacknifed 
MLE estimates. 
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The variance reduction, or lack of it, for each 
estimator when it is jacknifed can be determined by comparing 



the computed 


variances 


that result 


from 


a 


given 


subsample 


size and 


the 


variance 


coefficients 


» » 


a 2 


, etc . 


If, for 


example , 


the 


variance 


coefficients 


( a q > 


a 2 


» a 3 » 


and 014 ) for 



each of the estimators (Figures 4a through 4d) are used, to 
determine the variance for a subsample sice of 500, the vari- 
ances for each estimator would be: 



VAR (MOMENT) = 


58 . 37 


174.72 

+ 


1500.13 


5446.25 

+ 


= .1273 


500 


500 1 ' 5 


5 00 2 


5 00 2 ‘ 5 


VAR (MLE) 


67.52 


961.91 


10183.5 

+ 


27503.3 


= .0848 




500 


500 1 ' 5 


500 2 


5 00 2 ' 5 




VAR (JNK MNT ) = 


79.20 


457.44 


5216.3 

+ 


14533.1 


= .1357 


500 


500 1 ' 5 


5 00“ 


5 00 2 * 5 


VAR (JNK MLE) = 


62.66 


825.2 


+ 8879.79 


24484.2 


= .0827 


500 


5 00 1 ' 5 


500 2 


2 5 
500^ ' 3 



These computations show a lack of change in the variance 
when the estimators are jacknifed and agree with the statis- 
tical and graphical conclusions drawn earlier. 

The bias in the MLE and moment estimators is difficult 
to compare graphically (Figures 4a § 4b) since the unbiased 
estimate approaches the estimated expected value so quickly 
and because of the compacting caused by the range of the 
vertical scale. For this reason reduced graphics are produced 
(Figures 5a through 5d) , which the reader will recall 
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suppresses the outliers. Since the scale is reduced from a 
range of 15.28 (14.00 - (-1.28)) to 8.10 (9.14 - 1.04) or 
almost 50%, the evolution of the bias as subsample size 
increases becomes clearer, as Figures 5a § 5b demonstrate in 
comparison to Figures 4a § 4b. These graphs indicate that, 
for at least small subsample sizes, the moment estimator is 
a more biased estimator than the MLE estimator. 

Reduced graphics are not of much use for the jacknifed 
estimators (Figures 5c 5 5d) since the bias is almost com- 
pletely removed by jacknifing. They are useful only in the 
sense that they are on the same scale as their nonjacknifed 
forms and a comparison of bias between estimators can only 
be done graphically if the graphics of estimators are on the 
same scale. This, however, was already accomplished from 
Figures 4a through 4d. 

As was done for the variance, the amount of bias for 
a given subsample size can be estimated by dividing the coef- 
ficients, 8^, etc., by the sample size squared, cubed, 

etc. , as appropriate and summing (See Chapter II "USAGE” for 
the mathematical formula) . The computations will not be 
carried out here as in the case of the variance but it is 
useful to compare the 3^ coefficients of the estimators 
since these coefficients are likely to contribute the most 
significant amount of bias (they are divided by the smallest 
denominator). In the case of the moment estimator, 3]_ is 
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16.50 and for the MLE estimator it is 14.77. This suspected 
greater bias by the moment estimator can be seen visually 
from the graphics. The smaller coefficients for 3^ of the 
jacknifed moment and MLE estimators are 1.47 and 4.09 respec- 
tively. This gives an indication of how much less the bias 
becomes when the estimators are jacknifed. This kind of 
comparison of 3 ^ coefficients should be done with caution, 
though, since subsequent coefficients C 3 2 » £ 3 * ••) may be so 
large that they contribute more to the bias than the 3 -^ coef- 
ficient. It should also be pointed out that these 3-^ coeffi- 
cients should be zero because, as mentioned earlier, the 
jacknife removes the 3^(l/n) term. In fact, their estimated 
standard deviations are 4.00 and 5.67 respectively, so that 
they are within one standard deviation of zero. 

As important as the size of these bias coefficients 
are their accuracy. Comparing the bias coefficients from 
Figures 5a through 5d in relation to their standard deviation 
with those of Figures 3a through 3d (Table 4) which were 
determined from only two replications , shows the advantage 
of large replication in estimating the bias of an estimator. 
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Table 5: Bias coefficients for each estimator's 

population of estimates and their standard 
deviations for small and large replications 
(M=2 and M-100) 





BIAS 

COEFF. 


NUMBER 


OF REPLICATIONS 


ESTIMATOR 


TWO (2) 


ONE HUNDRED (100) 




B 1 


7.40 

(13.5) 


16.49 
(2 .59) 


MOMENT 


B 2 


183.1 

(1105.5) 


-182.2 

(204.1) 




B 3 


1947.9 

(21284.5) 


3659.04 

(4526.5) 




*1 


32.65 

(19.1) 


14.77 

(2.83) 


MLE 


b 2 


-1998.7 

(1031.2) 


-14.12 

(220.2) 




B 3 


47308.0 

(18073.1) 


1602.74 
(4798 . 3) 




B 1 


6.59 

(11.0) 


1.47 

(4.01) 


JNK MNT 


B 2 


-912 . 54 
(1016.4) 


9.75 

(335.5) 




B 3 


21772 . 7 
(23318.6) 


-1369.26 

(7411.7) 




B 1 


75 . 33 
(16.23) 


4.09 

(5.67) 


JNK MLE 


B 2 


-5815.6 

(1450.6) 


-203.3 

(436.02) 




B 3 


112632.0 

(32037.1) 


793.11 

(9243.2) 
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SAMPLE SIZE (Nls 500 NO. Of REPLICATIONS (M»: 100 DECREE Of REGRESSION (Oil 
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ESTIMATORS MOMENT ESTIMATOR (RECIPROCAL Of SQUARED COEFFICIENT Of VARIATION! Of THE SHAPE PARAMETER OF THE GAMMA D 1 STR UUTi ON, 
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ESTIMATOR; MAXIMUM LIKELIHOOD ESTIMATE OF THE SHAPE PARAMETER OF THE GAMMA DISTRIBUTION 



SAMPLE SIZE INIs 500 NU. OF REPLICATIONS <M I ; 100 DEGREE OF REGRESSION ID): 



CO 

o 



i i 

►* i — * 

i i 



1 




om'O 


UJUJ 




1 




| (N-* 


<-•— *ao 


.in 


1 




UJ • • 


»*f*»o* 


u 


1 




v*<NinoO 


3>com 


0 * 


1 


0 


CO'OO* 1 


>lAO 




1 0 


0 


0*<n>0 


O'^'* 


• 


— ♦ * — *n 


— ♦ 


•com 


•— m 


0 * 


1 <N 




* • • 


* • • 





♦***oo« 



► ♦ ♦ ♦♦****ooo+- 



•+ooo* 



■*oot** 



OUTGO 

<Nm® 

— • mm *m • 

Of^'f 

I O>0 O'*-* 



”♦0000 
500 I 
J*'*'* 

•mm 

• ♦ 



• • min'* 





OOOOO 


>Q *<N 


-0 


O 


33 or- 1 


'OtnN 


CO 


•n 


a* -ro 


r-f»m 


• 


<N 


•mm 


'*m • 


*0 




* • t 


r-cgtfi 

•— *m 
o*-*m 


X 

N 

m 



-•Oh* 

00 m 

I <N.x 

uj • • 
lftN0»OO 
CT'vO'* I 

•-om 

* • • 



Ovom 
1 m-* 

hi • • 

»«<K\JOO 

m oiM'O 

fM 3MAN 



i/\h* m co 
■OuA't m 
mr»o* >* 
<000 • 

-*som ** *0 

•<o • 



* • « 



-•vOO* 
O«'0 
1 <Nm 




UJ 

-J 

X 

2 : 

<x 

i/>UJ 



** 

I 



OO 

I I 



mr-moo 


*• 




o<nco 






OO CD 






• *m 


Z 




m— « • 


UJ 

O 






X •• 


•• 




X</) •• 


lO 


-*®m 


LUXm 


X 


o**« 


321- 


2 


1 'T 'O 


UUJZ 


UJ 


UJ • • 


— UJ 


— * 


aom>*00 


1 u — 


O 


CT'eTVnO 


—u 


•• 


<7* (NO* 


c/TU. *— 


X 


• •m 


UllUl 


X 


'*-* • 


'JUJU. 


UJ 




<CUJ 


a 




acv JO 

UJ 0 


0 




> 1 
< 1 


1 


— .om 


z 


UJ 


omm 


20 2 


0 


1 **<N 


a-*a 


2 


UJ • * 


LO~ 


< 


(NnO®0-« 


2viw 


— 


a*r — 0 


Oujcn 


CC 




— acuj 


< 


• •** 


LOOOC 


> 


'**■• • 


LOUJO 






LUQCLU 


X 




cc ac 

OX 


0 




uicu. 


z 




X 3 


a 


2i/)tn 


UJ 




<tn— 


XU> 


in 


LULU OO 


OZuj 


LO 


120 


<xo 


UJ 


Z 3C>- 


2— 


a: 


<oauja 


<acQ 


0 


LUH-»*- XO 


jj<x 


UJ 




2><n 


X 



IN'* 
-*cd 
cd a* 
(NO* 
• • 
-*m 



11 n 

2X 

— < 



< 

o 



< 

u 



55 



ESTIMATORS FOUR FOLD JACKN1FED MOMENT ESTIMATOR OF THE SHAPE PARAMETER OF THE GAMMA DISTRIBUTION 

Figure 4c 
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SAMPLE SIZE (N|s 500 NU • OF REPLICATIONS IMI s 100 OEGREE OF REGRESSION (D|; 





O 




a 




+ 




— 




Ui 




Ui 




<r<o a 


A 


O 




O®"^ 


M 


cO 


-nrvr- 


• co * 


• 


mm 


OmO 






X 


l or- 


iflON 


<r 


t— 


Ui • • 


-ONX\ 




A 


rsir^f— 00 


m w 


A 


— 


| | 








O'*"-* 








• mm 






< 


A • • 






x 





UJ • • 
—«<M -TOO 


m^OvO 


m 


a 


0 


Nf CO<J» I 


mm *<M 


mm 




A 


OOA 


CMr**-+ 


• 


X 


<M 


• Arr> 


•«o • 


0 


Ui 




A • • 


(MvO<f* 


O 


»— 






ao-uo 


A 


UJ 






•UsTfM 


mm 


X 






1 


1 


< 










X 










< 




mrnrtym 

OO'O 






X 




1 m~* 






UJ 




UJ • • 






X 




-«7*NOO 






< 


-0 


ao*m 1 






2 


>0 


o<rr*- 






A 




• 'Om 

A • • 






UJ 






>COsT 


>0 


X 

>— 






AOO* 


H 








O'r-A 


r» 


X 






fUJCO 


• 


a 




~»0<7* 


• mD A 








OOH 


O • • 


r» 






1 -T-* 


-+0<M 


mm 


2 




■UI • • 






a 




®<M-*<PO 






— • 


A 


j*>rfM 






tap 


CM 


OfMsO 






< 


— 


• +-cn 










A * • 






X 










< 










> 






<M-* 




X 






OO 




0 




*^C0 

O'O'-* 


1 1 
UIUJ 








1 m<M 


om® 


m 


2 




UJ • • 


.Mm a 


A 


Ui 




f'COuOO 


OCO<M 


O 




O 


NO« 1 


®r-m 


m 


U 


O 


»-«OA 




• 




mmt 


• com 


•— «m 


® 


X 




A • • 


>r • • 


A 


X 



a 

u 



>qoooo+- 



b I- 
O 



7* 

o 



*1 • 

o 



- I - 

o 



•I * 

o 



>0 

o 



» I- 

o 



• I - 
— 
o 
+ 
Ui 

o 

o 



0-H>0 
I Am 
UI • • 



► I ■ 

o 



r-co— OO 


•• 




■Omrnmrn 


A 




— or- 






•a»m 


2 




A • • 


Ui 

O 






X •• 


•• 




U A *♦ 


A 


— 0*0 


UJh-A 


X 


OnT<M 


CZt" 


2 


1 >0® 


UUJZ 


Ui 


UJ • • 


— UJ 


— 


aa*too 


1 U-* 


O 


or—* 


-0 


— 


<M^I*» 


AVL** 


X 


• .m 


UiU.X 


X 


A** • 


OUJX 


X 




<QLU 


0 




xoo 

UI O 


0 




> 1 
<■ 1 


1 


-*AA 


2 


X 


OOA 


zcz 


0 


t O'® 


0-0 


2 


tu • » 


A— 


< 


— OCOO'- 


2AA 


— 


-•cocn 


OUJA 


X 


>r>rco 


— qcUJ 


-<1 


• • m 


iflOtf 


> 


A mm • 


AUJO 






Uiotuj 


z 




X X 
Ox 


0 




UiQX 


2 




x a 


a 


2</n/> 


X 


— 



x < a— 

3! UJUJA 

< £ZQ 

Loui 2 

23 fU <OOLUa 
O— UiX*-xO 
aa xaaa*: 



U-0> 

QZUi 

<0 

2— 

<XO 

Ui<l- 

2>A 



at 

a 



o 

at 

x 

>rm — 

CT(M O 

crv*r ui 

O— X 



— o 



X 

o 



HI) 2 

ZX I— 
vi 

XX UI 



2 

O 

X 



X 

o 



57 



Figure 5a (REDUCED GRAPHICS) 
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Figure 5b (REDUCED GRAPHICS) 
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Figure 5c (REDUCED GRAPHICS) 
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3 . Case 1, Example 3: Large Subsample Sizes and Large 

Replications (K=5.0) 

Figures 4a through 4d showed the necessity of increas- 
ing subsample sizes in order to examine the normality of the 
estimates. Figures 6a through 6d are for increased subsample 
sizes (100 to 2000) with a large number of replications (M=25). 
The box plots in these figures show the increase in normality 
with increasing subsample size n^ when they are compared with 
small subsample sizes (Figures 4a through 4d) . The reduction 
in variability is also obvious by noticing the lack of outliers 
in the box plots and the smaller range of the vertical scale 
for larger sample sizes. The standard deviation for each of 
the estimators for small subsample sizes (ng=500, Figures 4a 
through 4d) is compared to those for large subsample sizes 
(n^=2000, Figures 6a through 6d) , Table 6. It can be seen 
that the variance is reduced for larger subsample sizes by 
the appropriate factor of two. However, the degrees of freedom 
in the standard deviation of the estimators are 25 at n^ = 2000 
and 100 at n^ = 500. Again it is clear from these figures 
that jacknifing causes no lack of precision in either the MLE 
or moment estimators. 

The bias in the moment and MLE estimators is difficult 
to compare from the graphics because of the range of the verti- 
cal scale. For this reason reduced graphics (Figures 7a § 7b) 
for these two estimators are produced. Reduced graphics for 
the jacknifed estimators would not be very useful since it 
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has been shown in previous examples that very little bias 
exists. For this reason they are omitted. 

Table 6: Standard deviations of the population of 

estimates of the shape parameter for each 
estimator for subsample sizes 500 and 2000 



SAMPLE SIZES 



ESTIMATOR 


ng = 500 


n 6 " 20 


MOMENT 


. 3417 


.1592 


MLE 


. 2730 


.0910 


JACKNIFED 

MOMENT 


.3400 


.1597 


JACKNIFED 

MLE 


. 2759 


.1090 



A change seems to have occurred if a comparison is 
made of Figures 7a 5 7b with those of Figures 5a § 5b. The 
graphics seem to indicate that the MLE ' s unbiased estimates 
have shifted from being biased on the high side to the low 
side for subsample size n^=100 and possibly n 7 =200. Also, 
it appears that the two estimators have about the same bias 
or that possibly the moment estimator has less bias. Looking 
further at the bias coefficients that determine these graphics, 
it can be seen that the coefficients in Figures 5a and 5b are 
estimated with more precision than those of Figures 7a and 
7b. One may be tempted to assume that the MLE estimator is 
more biased for larger sample sizes. But this is not likely 
and the change in the bias structure for Figure 7b can be 
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attributed to the difference in the values of the bias 
coefficients (3-^, 3-?, ..) compared to those of Figure 5b. 

These coefficients should be the same but the coefficients 
of Figure 5b were computed from 100 estimates (100 regressions) 
and those from Figure 7b only 25. The precision can be 
expected to be better in the first case. This is obviously 
the case if the standard deviation of the coefficients are 
compared, and more confidence can be placed in the bias 
structure of Figure 5b. 
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SAMPLE SUE (Nls 2000 NO • Of REPLICATIONS (MJs 25 DEGREE Of REGRESSION 4DJ 
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sample SIZE (NIi 2000 NO. OF REPLICATIONS <M I : 25 OEGREE Of REGRESSION <0)2 
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ESTIMATORS MAXIMUM LIKELIHOOD ESTIMATE Of THE SHAPE PARAMETER Of THE GAMMA 01 STR itiUTION 
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VERTICAL SCALES YMIN = 2.5001 

YMAx = 7.7173 

ESTIMATORS FCUR FOLD JACKNI FED MOMENT ESTIMATOR OF THE SHAPE PARAMETER OF THE GAMMA DISTRIBUTION 



SAMPLE SIZE INI J 2000 NQ • UP REPLICATIONS INI* 2S DEGREE OP REGRESSION (D)« 
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ESTIMATOR^ FCUR POLO JACKNIPED MAXIMUM LlKELlHUOO ESTIMATE OP THE SHAPE PARAMETER OP THE GAMMA DISTRIBUTION 

Figure 6d 



SAMPLE SIZE (N): 2000 NLi • Of REPLICATIONS tMi : 25 OEGREE Of REGRESSION IOJS 
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SAMPLE SIZE (N|: 2000 NO. Clf* HEP L 1CA T IONS (Ml J 25 OEGREE Of KtuKfcSSlUN ( 0 1 • 



O' 

a 

+ 

x 

v^mco ♦ 

• O • -fl 

m mm 'O 
f>»axn -♦ 

a>wro 
-4 «~4 -• 

I 



o 

-o 

o 



1 

1 




N» 






1 




<3 






1 




♦ 






l 




UJ 






1 




'T'Of- 






1 




ma>>t 


« 




1 




• o • 


m 




1 




rx*(*v» * 


O' 




1 




— *n 


o 




1 




oaxo 


o 




1 

1 




-r *<m 
l 


i 




1 

1 

1 










1 


OOC'O 








1 


I l cun 








1 


XX • • 








1 








• 


1 


O'fON | 






2 


1 


OOwOJ 






O 


1 


O • O'w 






w 




rs»U> • • 






w 


1 








o 


1 




«*tna* 




33 


1 




>00 


in 


w 


1 




O— nj 


• 


aC 


1 




cm *a» 


-4 


w 


1 




• cm . 


— * 


00 


1 


omo 


33 « 


-4 


■4 


1 


i ■ *■ < 


rO» •» * 


m 


a 


1 


UJ « • 








1 


OmcDOO 






< 


1 


O'O'O'O 






x 


1 


OO vfO 






X 


1 


O *wf\| 






< 


1 


—•m ♦ • 






o 


1 

1 








UJ 


1 










1 




row 




►— 


1 




oo 






1 




1 1 




u. 


1 


o«o<* 


UJUJ 




a 


1 


1 c*><M 


ONO 


o 




1 


UJ • • 


CMOCT* 


OJ 


at 


1 


nooOO 


romO 


>r 


UJ 


1 o 




moo 




V- 


CM — O 


OONN 


O aor* 


• 


UJ 


1 o 


m •nj<M 


• evj.— < 


CO 


X 


1 «* 


in • • 


m • • 


rn 





1 


"•OO 










1 


omoo 








< 


I 


i xf'-r 








21 


1 


UJ • • 








to 


1 


J-'fl'TOO 


•• 








1 


ONmm 


oo 






UJ 


} 


m^oco 


b» 






X 


1 


cm *vrru 


2 






X 


1 


in • • 


UJ 








1 




■4 






X 


1 




O 






a 


1 












1 




u. •• 


H 




X 


I 




U.l/3** 






w 


1 


wwO 


X*-</» 


>— 




< 


1 




az)- 


2 




X 


1 


i -rnj 


UUJZ 


UJ 




•4 


1 o 


UJ • • 


WJJ 


■4 




•— 




>r<Mxroo 


1 U*- 


u 


0*0 


to 


1 m 


O'TV'O 


— u 


•4 


n in 


X 


1 


Owner* 


1/>U»— 


X 


CMro 




1 


CM ‘vTCM 


Luu.u„ 


X 


• • 


Q 


I 


in • • 


OUjU. 


Ul 


or- 


a 


1 o 




<aui 


◦ 




o 


— «n 




auO 


o 




X 


1 CM 




UJ <_) 






■4 


1 




> i 


1 




— r 


1 




< 1 




H H 


X 



l n*0 
x • ♦ 
O'OCOOO 

0-*rno 

0<M >0<7* 



2 

202 
C— O 

t/»— 

2</i vi 
Ox«o 
-*<xx 



2x 

— <f 
XX 



O 

x 



0*0 








♦ 








o ow 


1 o 

— o 

1 

1 

1 


-» •'OCVJ 

in • • 


too at 

lOLUO 

xauj 

esc a. 

OX 

XGX 

rf q 


> 

2 

O 

2 

Q 


X 

X 

<t 

o 

to 


X 

< 

X 

•• 

gf 




1 








• 








J. j. 


1 


m*4 


— 1 

w 




1 *■■4 
W 




w 


~ i. 




x 2toto 


X 




X 


O 


o o 


o 


o 


o 


o 


o 


o 


O 


o 




a. <t04- 


xo> 


to 


< 


X 


4 * 


♦ 


* 


♦ 


♦ 


♦ 


♦ 


♦ 


4> 




T. XX to 


0 2X1 


to 


<-> 


< 


X x 


X 


X 


X 


X 


X 


X 


X 


X 




< X20 


<o 


X 


4. 


X 


<r o' 


n 


4 


r*» 


m 


ao 


vT 


o 


o 




<OXZ 3W- 


2w 


X. 


V- 


"4 


P» vO 


o 


o 


n 


n 


vT 




-T 


ro 




DM<OQlLOt: 


<QCO 


o 


at 


►- 


• • 


• 


• 


• 


• 


• 


• 


• 


• 




O— UJW*— 


X<Il— 


X 


X 


to 


0 

0 


o 


o 


o 


o 


O 


O 


o 


o 




otoX totoio*: 


X>to 


X 


> 


X 



69 



Figure 7b (REDUCED GRAPHICS) 



B. EXAMINING THE GAMMA DISTRIBUTION FOR THE SHAPE PARAMETER 
(K=0 . 2 ) 

1 . Case 2, Example 1: Small Subsample Sizes and 

Small Replications (K=0.2)' 

In the preceding case for K=5.0 four estimators of 
the shape parameter of the gamma distribution were evaluated 
for varying subsample sizes and replications. This evaluation 
is repeated here with the shape parameter changed from K=5.0 
to K=0.2. This second set of comparisons is not only useful 
in demonstrating how the RAGE package shows the behavior of 
different estimators in regards to bias and variance, but 
also how these estimators may differ when the parameter being 
evaluated takes on different values . 

Figures 8a through 8d are graphical results of the 
four estimators of the shape parameter, K=0.2. For these 
small subsample sizes and small population of estimates it 
is apparent from the graphics that the MLE estimator (Figure 
8b) has tighter box plots than does the moment estimator 
(Figure 8a). This shows graphically that the variance for 
the MLE estimator is less than that of the moment estimator. 
This was not obvious when the shape parameter was 5.0. In 
Figures 3a § 3b of Case 1, Example 1, it was difficult to 
distinguish any difference in the box plots. This comparison 
shows how the MLE estimator is much better at reducing the 
variance when the actual shape parameter is as small as 0.2. 
Or, alternatively, it may be an indication of how much xvorse 
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the moment estimator is for estimating shape parameters 
that are small. 

Figures 8c 8 8d are the jacknifed forms of the esti- 
mators. It is not at all clear from the graphs that the 
jacknife is affecting the variance of the estimates. This 
lack of change could be attributed to the small number of 
estimates which suggest that the number of replications 
should be increased to achieve a better graphical evaluation. 

In Case 1, Examples 1, 2 and 3, the coefficients of 
variance (a-p and a^) were helpful in determining 

which estimators produced the least variance. A look at 
in Figures 8b and 8d shows that they are negative numbers. 

Since often contributes the most significant amount to 
the variance, there is some suspicion that a determination 
of the variance by the use of these coefficients could result 
in negative variances for the MLE estimator and its jacknifed 
form. If a ? and are not large positive numbers (relatively) 
it is especially likely. In these two cases a 7 is a rela- 
tively large positive number and negative variances do not 
occur if the calculations are carried through. However, 
for lesser degrees of regression (i.e. D=2) negative variances 
can occur. This is a result of the wide fluctuation of the 
variances of the small subsample sizes whose values are used 
as dependent variables in the regression which determine 
otp a 9 , and ou . These variance coefficients are often 
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questionable when the variances of the population of esti- 
mates vary over a wide range and further suggest increasing 
the number of replications so that the variance coefficients 
can be more accurately determined. 

The bias in the moment estimator (Figure 8a) appears 
to be considerably more than that of the MLE estimator 
(Figure 8b) when the graphs of figures are compared. This 
bias is further reduced when the jacknife is applied to the 
estimators (Figures 8c 8 8d) . A comparison of the bias is 
inadvisable at this point since the bias coefficients have 
large standard deviation values in comparison to their 
estimates. The 3^ coefficient, for example, of the MLE 
estimator is .0864 and its standard deviation is .189. As 
mentioned in previous examples, 3^ will normally contribute 
the most to the bias and since its standard deviation is 
larger than the estimate itself, its accuracy is highly 
questionable. Before any firm conclusions can be drawn 
concerning the bias of the estimators more estimates are 
required, which implies increasing the number of replications. 
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ESTIMATORS MOMENT ESTIMATOR IRECIPROCAL Of SQUARED COEffICItNT Of VARIATION! Of THE SHAPE PARAMETER Of THE GAMMA OISTRlbOTION 

Figure 8a 
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2 . Case 2, Example 2: Small Subsample Sizes and 

Large Replications (fc=0.ZT 

Figures 9a through 9d are for small subsample sizes 
but large replications (M=100) . This increased number of 
estimates produces more definitive graphics for better 
evaluations of the bias and variance of the estimates pro- 
duced by the estimators. 

The graphs of the moment and MLE estimators (Figures 
9a § 9b) show clearly as before that the MLE estimator has 
much less variance than does the moment estimator. Both, 
however, have positive skewness and high nonnormal estimates 
up to subsample size n_j = 125. The jacknifed form of the 
moment estimator has about the same variance as its nonjack- 
nifed form but is more symmetrically distributed. This can 
be seen in a comparison of the box plots of Figures 9a § 9c. 

The distribution of the population of estimates for 
the jacknifed MLE estimator shows little if any difference 
graphically or statistically than the nonjacknifed estimates. 
This can be seen by comparing Figures 9b S 9d , which show 
little difference when the box plots are compared for 
various subsample sizes. For large values of the shape 
parameter, i.e. K=5.0, it could be seen from Figures 4b § 4d 
of Case 1, Example 1 that the jacknifed MLE estimator 
increased the symmetry of the estimates. This was obvious 
both in the graphics and by comparing the skewness coeffi- 
cients. For K=0.2 the only indication that the estimates 
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are more symmetrically distributed is the comparison of the 
skewness coefficients for small subsample sizes (n^=33, 
n2=50 and nj=71) of the MLE estimator and its jacknifed form 
(Figures 9b 5 9d) . The skewness appears to decrease slightly, 
but not discernably, for larger subsample sizes. 

The amount of variance that was observed graphically 
can also be seen statistically by computing the variance of 
each estimator for designated subsample sizes using the vari- 
ance coefficients , a~ and of each estimator. These 

coefficients for each estimator will result in positive 
variances. This is not always the case when the coefficients 
are determined from subsample sizes with a small population 
of estimates. This example illustrates the increased confi- 
dence that can be put into the variance coefficients as the 
number of estimates (replications) increases. 

The bias again appears to be substantially greater 
for the moment estimator than the MLE estimator for comparable 
subsample sizes. Here, because of the range of the vertical 
scale, reduced graphics are produced for the moment and MLE 
estimators to better compare the bias of these estimators 
by reducing the scale (Figures 10a and 10b) . A comparison 
of Figures 10a and 10b clearly shows the increased bias 
caused by the moment estimator. It was shown in Example 1 
of Case 1 that the jacknifed estimators have little bias 
remaining in the estimates . This appears to be the case for 
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K=0.2, as indicated by the small bias coefficients of the 
jacknifed estimators (Figures 9c § 9d) . For this reason, 
reduced graphics for the jacknifed estimators are omitted, 
since indications are they would show little more graphically 
than Figures 9c § 9d. 

It is interesting to note that the bias of the moment 
estimator compared to that of the MLE estimator is not only 
greater, but greater for small values of the shape parameter 
(i.e. K=0.2) than for large values (i.e. K=5.0). This can 
be seen by comparing the bias in Figures 10a 5 10b with that 
of Figures 5a 5 5b. This indicates that the moment estimator 
is not only more biased than the MLE estimator but consider- 
ably more biased when the shape parameter being estimated is 
small. 
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ESTIMATOR; MOMENT ESTIMATOR IRECIPRQCAL OF SQUARED COEFFICIENT OF VARIATION! OF THE SHAPE PARAMETER OF THE GAMMA DISTRIBUTION 
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ESTIMATOR: MAXIMUM LIKELIHOOD ESTIMATE Of THE SHARE PARAMETER OF THE GAMMA DISTRIBUTION 

Figure 9b 
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ESTIMATOR: FOUR FOLD JACKNIFEO MOMENT ESTIMATOR OF THE SHAPE PARAMETER OF THE GAMMA DISTRIBUTION 

Figure 9c 
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ESTIMATOR; FCUR FOLD JACKNIFEO MAXIMUM LIKELIHOOD ESTIMATE OF THE SHAPE PARAMETER OF THE GAMMA DISTRIBUTION 

Figure 9d 
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ESTIMATORS MCMENT ESTIMATOR (RECIPROCAL OF SQUARED COEFFICIENT OF VARIATION! OF THE SHAPE PARAMETER OF THE GAMMA DISTRIBUTION, 

Figure 10a (REDUCED GRAPHICS) 
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ESTIMATORS MAXIMUM LIKELIHOOD ESTIMATE Of THE SHARE PARAMETER Of THE GAMMA DISTRIBUTION. 

Figure 10b (REDUCED GRAPHICS) 



Case 2, Example 5: Large Subsample Sizes and 

Large Replications (K=oT2)~ ~ 



o . 



As a last example, graphs for large subsample sizes 
with a large number of replications (M=25) are produced for 
each estimator. The box plots of these graphs (Figures 11a 
through lid) look similar to their counterparts for small 
sample sizes and large replications (Figures 9a through 9d) . 
This is due to the difference in the vertical scale and does 
not indicate a lack of variance reduction by larger sample 
sizes. The range of the vertical scale of Figures 11a 
through lid is .495 and is smaller than that of Figures 9a 
through 9d , which is .9865, by a factor or two. This is a 
simple but important illustration of the caution that must 
be exercised when comparing graphics that are not of the 
same scale. 

Other than the reduced variance caused by larger 
sample sizes, the same conclusions that were reached in 
Example 2 of this case concerning the variance and the bias 
of the estimators is shown to hold here. The MLE estimator 
shows considerably less variance and bias than the moment 
estimator. The jacknifed moment estimator results in more 
normally distributed estimates and less bias. The jacknifed 
MLE estimator shows some improvement toward normality of its 
estimates and its variance has changed little. 

As in Case 1 for K=S.O, the larger the subsample 
size the more normal are the population of estimates. This 
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increased normality can be seen graphically by comparing the 
graphs of Figures 11a through lid with those of Figures 9a 
through 9d. The lack of outliers on the high side of the 
box plots for Figures 9a to 9d is an indication of convergence 
normality . 

Table 7 summarizes the results of the comparisons of 
the four different estimators for K=5.0 and K=0.2. A com- 
parison of the standard deviations shows the moment estimator 
to be consistently higher than the MLE estimator. These 
standard deviations, however , show little change when the 
primary estimators are jacknifed, indicating little or no 
change in variance by jacknifing. The bias values of the 
moment estimator are consistently higher and the root mean 
square error (r.m.s.e.) is considerably larger for K=0.2 in 
comparison to the MLE estimator, indicating how poor the 
moment estimator is for evaluating shape parameters of gamma 
distribution that have small values. 

These conclusions were easily drawn from the graphs 
and the statistical output provided by the RAGE program and 
did not require extensive programming and analysis by the 
experimenter. This kind of expedient examination of estima- 
tors not only saves the experimenter time but reduces signifi- 
cantly the computer time required for analysis. This 
facilitation of analysis is the primary purpose of the RAGE 
methodology . 
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Table 7: Bias and variance comparisons for varying 

subsample sizes and shape parameters 

SAMPLE § 

SHAPE SUBSAMPLE 



PRMTR 


SIZE 


M 


ESTIMATOR 


ESTIMATE 


STD. DEV. 


BIAS R 


. M . S . E 


1 

Ol 

o 


500 


100 


MOMENT 


. 2071 


. 002927 


-.0071 


.03012 








MLE 


.2006 


. 0008701 


- .006 


. 01057 








JN MOM 


.2015 


.003211 


- . 0015 


. 03215 








JN MLE 


.2000 


. 008707 


.00 


.0087 


0.2 


2000 


25 


MOMENT 


.2021 


.003584 


-.0021 


.01804 








MLE 


.2001 


.0008417 


.0001 


.00421 








JN MOM 


.2005 


.003644 


- .0005 


. 01823 








JN MLE 


.1997 


.0008966 


. 0003 


.00494 


5.0 


500 


100 


MOMENT 


5.012 


.03417 


- .012 


.3419 








MLE 


5.095 


. 02730 


- . 095 


.2891 








JN MOM 


4.981 


.0340 


.019 


. 3405 








JN MLE 


5.057 


. 02759 


- .057 


.2817 


5.0 


2000 


25 


MOMENT 


4.994 


. 03184 


.006 


.1593 








MLE 


5.049 


.01821 


-.049 


.1033 








JN MOM 


4.988 


.03194 


.012 


.16015 








JN MLE 


5.014 


.02180 


- .014 


.1099 
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VERTICAL SCALE; YMIN * -0.0744 

YMAX «= 0.4206 

ESTIMATORS MOMENT ESTIMATOR (RECIPROCAL OP SOUAREU COEPPICIENT OP VARIATIONI Op THE SHAPE PARAMETER OP THE GAMMA 0 1 STR I bUT I CJN 

Figure 11 a 
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VERTICAL SCALE: VMLN ; -g.gjft 

ESTIMATOR: MAXIMUM LIKELIHOOD ESTIMATE OP THE SHAPE PARAMETER OP THE GAMMA DISTRIbUTION 

Figure lib 
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ESTIMATOR: FOUR FOLU JACKN1FEC MOMENT EST IHAIOR Of THE SHAPE PARAMETER OF THE GAMMA DISTRIBUTION 

Figure 11c 



SAMPLE SIZE (N): 2000 NU • UP REPLICATIONS (MI; 25 IfEGREE OP REGRESS |UN IU); 
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